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SPECTRAL  ANALYSIS  VIA  LAG-RESHAPING  OF  THE  CORRELATION  ESTIMATE; 
PROGRAMS  AND  SIMULATION  RESULTS 


INTRODUCTION 


The  fundamental  performance  of  the  generalized  spectral  analysis 
technique  employing  quadratic  frequency-smoothing  of  Fourier-transformed, 
overlapped,  weighted  data  segments  was  thoroughly  investigated  analytically 
and  reported  in  [1].  One  attractive  possibility  pointed  out  there  was  that  of 
doing  lag-reshaping  of  the  first-stage  correlation  estimate,  and  thereby 
realizing  effective  spectral  windows  with  low  sidelobes  and  good  decay  rates, 
without  the  need  for  overlap  or  any  temporal  weighting  at  all.  However,  no 
programs  or  simulation  results  for  this  particular  technique  were  presented  at 
that  time.  We  rectify  this  situation  here  by  presenting  programs  for 
achieving  maximally  stable  low-sidelobe  spectral  estimates  via  the 
lag-reshaping  procedure  [1,  pages  36-40]  and  then  employ  these  routines  to 
exhibit  some  simulation  results  for  data  with  tones  and  noise. 

Spectral  analysis  techniques  have  received  a  great  deal  of  attention  in 
the  past' [2  -  13],  ranging  from  the  original  autocorrelation  approach  of 
Blackman-Tukey  [3]  to  the  more  recent  weighted,  overlapped,  segment-averaging 
FFT  approach  [8  -  13].  These  two  apparently  disparate  approaches  are  limiting 
special  cases  of  a  generalized  framework  [1]  for  spectral  analysis;  thus 
consideration  of  this  general  technique  elucidates  the  fundamental  behavior 
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and  performance  of  a  rather  wide  variety  of  spectral  approaches  and  their 
tradeoffs.  This  generalized  framework  was  originally  presented  in  [14  -  16], 
where  a  brief  summary  of  some  of  the  main  features  was  mentioned.  The 
analytical  results,  detailed  derivations,  and  quantitative  results  were  given 
in  [1]. 


There  are  two  fundamental  parameters  that  critically  affect  the 

performance  of  any  spectral  estimation  technique.  They  are  the  available 

record  length,  T,  of  the  stationary  random  process  under  investigation,  and 

the  effective  frequency  resolution,  Bg,  of  the  technique  under  consideration 

We  would  like  to  be  able  to  attain  fine  resolution  (small  B  )  with  short 

e' 

data  lengths  and  storage  (small  T);  however,  stable  results  (small 
fluctuations)  are  achievable  only  if  the  product  TBg  is  much  larger  than 
unity.  The  problems  of  interest  are  how  to  make  optimum  use  of  a  given 
limited  amount  of  data,  in  order  to  realize  a  specified  desired  resolution 
with  maximum  stability,  and  to  determine  what  tradeoffs  are  available 
regarding  windowing  and  weighting  at  different  stages  of  the  spectral  analysis 
procedure.  It  is  assumed  that  the  reader  is  familiar  with  the  tradeoffs 
presented  in  [10]  for  the  weighted,  overlapped,  segment-averaging  FFT 
procedure,  and  with  the  concepts  and  results  of  the  generalized  framework  in 
[!]• 


In  this  report,  we  will  confine  attention  solely  to  the  case  of  abutting 
rectangular  temporal  weightings  with  no  overlap;  this  procedure,  with 
appropriate  lag-reshaping,  has  been  shown  [1,  pages  48-50  and  58]  to  have 
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excellent  effective  spectral  windows  (low  sidelobes  and  rapid  decay)  and 
virtually  ideal  variance  reduction  capability,  under  proper  choice  of  the  lag 
weighting  and  the  ratio  of  weighting  lengths.  Furthermore,  it  requires  no 
temporal  multiplications  whatsoever  on  the  given  data,  and  it  avoids  the 
additional  manipulations  associated  with  overlapped  data  segments.  Thus  it  is 
a  strong  candidate  for  consideration  in  spectral  analysis. 

Separate  programs  for  spectral  analysis  of  complex  data  as  well  as  real 
data  are  presented  here.  Also,  the  programs  are  written  in  such  a  fashion 
that  if  multiple  frequency-resolution  bandwidths  are  desired,  they  can  be 
easily  accommodated  without  re-doing  the  bulk  of  the  required  data  processing; 
in  fact,  one  new  lag  weighting  and  one  FFT  will  suffice  for  each  different 
required  resolution. 

Although  the  present  programs  and  study  are  limited  to  auto-spectral 
analysis,  they  can  be  easily  generalized  to  incorporate  cross-spectral 
analysis  if  desired,  both  for  complex  as  well  as  real  data.  The  present 
conclusions  on  execution  time  and  stability  should  carry  over  directly  to  this 
more  general  case. 
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ULTIMATE  STABILITY  ATTAINABLE  FROM  A  GIVEN  RECORD  LENGTH 

Suppose  a  stationary  (complex)  data  record  x(t)  of  length  T  seconds  is 
available  and  that  we  wish  to  estimate  its  power  density  spectrum  G(f)  at 
frequency  f,  with  an  effective  frequency  resolution  of  Bg  Hz,  where  WQ( f ) 
is  the  narrowband  window  through  which  the  power  density  spectrum  is  to  be 
observed.  These  two  frequency-domain  quantities  are  related  according  to 


B 


(1) 


e 


This  bandwidth  measure,  B  .  is  called  the  statistical  bandwidth  of  W  (f) 
in  [6,  page  265].  The  relation  of  effective  bandwidth  Bg  to  the  half-power 
bandwidth  is  considered  in  (1,  appendix  A];  it  is  shown  that  for  good  windows, 
the  ratio  of  the  two  bandwidths  is  relatively  independent  of  the  exact  window 
shape.  Thus  it  is  possible  to  translate  results  to  other  bandwidth  measures 
without  significantly  affecting  the  essential  quantitative  aspects. 

If  we  take  the  original  data  record  and  pass  it  through  a  narrowband 
linear  (complex)  filter  with  power  transfer  function  equal  to  the  window, 
|H(f)p  =  WQ(f),  and  which  is  centered  at  a  frequency,  f  ,  of  interest, 
we  will  have  lost  no  relevant  information  about  the  process  in  the  frequency 
band  of  interest,  because  we  have  filtered  out  information  of  no  use.  We  can 
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now  estimate  the  power  in  the  narrowband  filter  output  process  and  use  it  as  a 
measure  of  the  spectrum  of  the  input  process  in  the  neighborhood  of  frequency 
f  ;  see  figure  1 . 


Figure  1.  Power  Transfer  Function  of  Narrowband  Linear  Filter  H(f) 


A  quality  ratio  can  be  defined  for  the  power  estimate  P  at  the  filter 
output  as 


(2) 


Under  the  assumptions  that  the 

filter  input  spectrum  G(f)  does  not  vary  rapidly  with  respect  to  bandwidth 

B  ,  that  the  observation-time  resolution-bandwidth  product  TB  is  large, 

6  ° 

and  that  the  filter  output  is  approximately  Gaussian,  it  is  shown  in 
[1 ,  pages  3-5]  that 

0  =  •  (3) 


n  _  Var(P)  P2  -  P 
U  -  2  - 

Av( P)  _2 

P 

where  the  overbar  denotes  an  ensemble  average. 
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This  is  the  smallest  value  of  quality  ratio  Q,  for  specified  resolution  B 

e 

and  available  record  length  T;  no  other  spectral  analysis  procedure  can 
outperform  this  benchmark.  Thus  (3)  serves  as  a  very  useful  comparison 
standard  for  any  technique  and  will  be  employed  here  in  the  stability 
investigation.  Specifically,  the  normalized  quality  ratio  (NQR)  for  a 
particular  technique  is  defined  as  the  quotient  of  the  quality  ratio  of  that 
technique,  (2),  relative  to  the  minimum  value  in  (3).  Thus  NQR  is  always 
greater  than  unity,  with  small  values  being  desirable. 
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SPECTRAL  ANALYSIS  PROCEDURE 


The  available  stationary  data  process  x(t)  is  presumed  to  be  sampled  at 
unit  time  increment,  for  convenience;  that  is,  =  1  second  in  [1,  (24)  et 
seq.].  There  should  be  no  problem  in  scaling  these  results  to  a  general  time 
sampling  increment  At.  Thus  the  observation  time  T  measures  both  the 
available  data  record  length  and  the  number  of  available  data  points. 

This  total  of  T  data  points  is  broken  into  K  abutting  non-overlapping 
data  segments  each  of  length  L1 ;  thus 

T  -  K  l1  ,  (4) 

where  L^  is  the  length  of  the  (rectangular)  temporal  weighting.  The  data  points 

in  the  k-th  piece,  1  <  k  <  K,  are  labeled  for  0  <  n  <  L]  -  1,  where  n 

is  a  time  index;  thus  in  terms  of  original  process  x(t),  we  have  (since  A^.  =  1 ) 

x^  =  x(n+(k-l )L|)  for  0  <  n  <  L1  -  1  ,  1  <  k  <  K  .  (5) 

Notice  there  are  no  common  data  points  in  any  of  the  different  pieces. 


FIRST-STAGE  ESTIMATES 

According  to  the  procedure  described  in  [1],  we  perform  a  forward 
N-j -point  (first-stage)  FFT  of  each  piece  of  data  (N]  =  power  of  2): 
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h-1 

Xmk^  =  2>  exp(-i2irnm/N1 )  x^  for  0  <  m  <  N1  -  1 ,  1  <  k  <  K,  (6) 

n=0 

where  m  is  a  frequency  bin  index.  Observe  that  no  temporal  weighting  is 
employed  on  the  available  data.  We  presume  that  FFT  size  N-j  and  segment 
length  are  chosen  to  satisfy  constraint 

^  >  h  •  (7) 

However,  we  will  discover  shortly  that  must  be  chosen  even  larger  than 

constraint  (7);  thus  some  zero-filling  is  required  in  (6),  namely, 

zeros,  prior  to  taking  the  FFT  dictated  by  (6).  The  spacing  of  the  frequency 

domain  components  {x^}  in  (6)  is 


and  they  cover  a  total  band  of  l/at  =  1  Hz.  Thus  the  m-th  frequency 
component  in  (6)  is  at  frequency  m/f^  Hz. 


(8) 


The  frequency  components  in  (6)  are  now  subjected  to  a  magnitude-square 
operation  and  segment-averaging  over  the  available  K  pieces  of  data,  yielding 
the  first-stage  spectral  estimates 
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(Scale  factors  will  be  accounted  for  later,  at  the  end  of  the  spectral 
procedure  description.)  These  power  spectral  estimates  also  occur  at  the 
frequency  spacing  indicated  in  (8),  and  cover  a  total  band  of  1  Hz. 

We  now  inverse  Fourier  transform  (9)  back  into  the  lag  domain,  obtaining 
first-stage  correlation  estimates  defined  as 

t1  .. 

an  =  >  exp(i2irnm/N1 )  Am  for  all  n.  (10) 

m=0 

This  is  an  N^-point  transform;  however,  we  let  it  define  the  sequence  {an"| 
for  all  n,  with  period  N1 .  The  spacing  of  correlation  estimates  ja  in 
the  lag  domain  Is  ■  1  second.  A  typical  representati ve  plot  of  {a^  in 
figure  2  reveals  an  important  property  of  correlation  estimates  obtained  via 
the  forward-and-inverse  discrete  Fourier  transform  procedure  of  ( 6) — ( 1 0) : 

GL  =  Sum 


Figure  2.  Wrap-Around  in  the  Lag  Domain 
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Namely,  {an}  equals  the  periodic  correlation  of  a  data  sequence  of  length 
L-j ;  unless  >  2L^ ,  {a^  will  suffer  wrap-around,  since  each  an  is 
the  sum  of  all  the  displaced  periodic  versions  of  the  desired  aperiodic 
correlation  (the  triangle  centered  at  n  =  0) .  However,  we  can  still  tolerate 
some  wrap-around,  if  the  lag  weighting  length  is  chosen  small  enough,  as 
indicated  in  figure  3: 


Figure  3.  Allowed  Lag  Weighting  to  Suppress  Wrap-Around 


Thus,  if 


L2  -  N1  '  L1  *  i,e“  N1  -  L1  +  L2  * 


(11) 


then  the  lag  weighting  wn  goes  to  zero  before 
from  the  aliasing  lobes  centered  at  n  =  ±N^ . 
correlation  estimates  via  an  FFT  approach  was 


any  undesired  overlap  occurs 
This  ability  to  get  unaliased 
also  pointed  out  and  utilized  in 
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[1 ,  page  12,  (41)  et  seq.]  and  [17 ,  (15)].  The  constraint  in  (11)  will  be 
developed  even  further  when  we  discuss  execution  time  and  stability  of  the 
final  spectral  estimate. 

The  constraint  in  (11)  is  not  an  upper  bound  on  and  l_2,  as  much  as 
it  is  a  lower  bound  on  N1 .  That  is,  if  we  want  to  estimate  the  correlation 
of  the  input  data  via  this  forward-and-inverse  FFT  approach  (for  minimum 
execution  time  reasons),  we  must  choose  the  first-stage  FFT  size,  N-j  >  +  l_2, 

in  order  to  circumvent  the  inherent  wrap-around  associated  with  the 
procedure.  By  contrast,  recognize  that  if  we  calculated  the  correlation  of 
the  available  input  data  by  the  brute  force  delay-multiply-add  procedure 
(Blackman-Tukey) ,  no  such  limitation  would  arise;  we  would  simply  compute 
correlation  estimates  to  the  maximum  lag  L of  interest.  Since  finer 
frequency  resolution  of  the  final  (second-stage)  spectral  estimates  is 
achieved  by  making  lag  L2  larger,  the  size  of  the  first-stage  FFT  will 
have  to  be  increased  accordingly;  in  fact,  it  must  accommodate  the  largest 
of  interest  in  the  second-stage  estimation  procedure. 

SECOND-STAGE  ESTIMATES 

The  second-stage  correlation  estimates  are  obtained  by  lag  weighting  the 
first-stage  estimates  in  (10): 

bn  =  wn  an  for  0  <  n  <  *-2  <  4  •  02) 
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Here  lag  weighting  coefficients  wp  are  selected  according  to 


w 


—  for  0  <  n  <  L2  <  L1  , 

17 


(13) 


n 


where  the  numerator  is  a  desirable  lag  weighting  with  low  sidelobes  and 
adequate  decay,  while  the  denominator  is  the  autocorrelation  of  the 
rectangular  temporal  weighting  employed  in  transform  (6);  see  [1,  (36)-(40)]. 

Actually,  since  transform  (10)  will  be  accomplished  via  an  FFT,  the 


stored  and  available,  while 


values  for  n  <  0  will  not  be.  Accordingly,  the  weighting  in  (12)  and  (13) 
must  be  augmented  to  account  for  a  symmetric  weighting  about  n  =  0.  This  is 
accomplished  by  also  weighting  the  upper  end  of  the  [a^  sequence  near 
n  =  N.| ,  as  indicated  by  the  dashed  curve  in  figure  3.  The  same  constraint 
(11)  suffices  to  guarantee  no  effect  due  to  aliasing  from  the  tail  of  the 
correlation  lobe  which  is  centered  at  n  =  0. 

Although  we  only  need  to  compute  {an}  up  through  L 2,  since  lag 
weighting  wn  is  zero  beyond  there,  the  FFT  will  yield  values  of  an 
via  (10).  The  excess  values  not  required  in  figure  3  are  simply  discarded. 
Also,  advantage  can  be  taken  of  the  conjugate  symmetry  of  an  about  n  =  0, 


±N, ,  etc.;  i.e.,  a  =  a„  since  a„  is  a  correlation  estimate, 
i  -n  n  n 

The  application  of  lag  weighting  (13)  in  (12)  can  be  accomplished  in  two 
steps: 
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a 

bn  =  wf,d)  n  for  0  <  n  <  L2  <  L,  .  (14) 

1 "  Li 

The  reason  for  this  separation  of  effort  is  that  the  division  of  ar(  by 

1  -  n/L1  can  be  done  once,  prior  to  knowledge  or  selection  of  l_2  and  w^  , 

stored,  and  then  re-used  for  several  different  desirable  lag  weightings  wj^  , 
which  can  have  different  lengths  L2  and/or  sidelobe  character  and  decay,  for 
whatever  frequency  resolution  is  of  interest.  Good  candidate  lag  weightings 
in  this  respect  are  the  optimum  ones  presented  in  [18].  The  restriction  that 
L2  <  L.|  in  (14)  is  necessary  to  avoid  a  division  by  zero  at  n  =  l2. 

Finally,  the  second-stage  spectral  estimates  are  obtained  by  transforming 
(12)  into  the  frequency  domain  (N2  =  power  of  2): 

L2 

Bm  =  exp(-i2irnm/N2)  bn  for  0  <  m  <  -  1  .  (15) 

n=0 


The  size  of  this  final  FFT  dictates  the  frequency  spacing  of  these  spectral 
estimates,  namely, 


L  =  L  H2 

M2 


N2At 


(16) 


Thus  the  m-th  frequency  component  in  (15)  is  at  frequency  f  =  m/N„  Hz. 

m  2 

In  order  to  be  able  to  observe  all  the  detail  of  the  second-stage  estimate,  we 
should  choose 

N2  »  L2  .  (17) 
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However,  N2  has  no  effect  on  the  stability  or  resolution  capability  of 
spectral  estimates  {b^;  rather,  the  stability  depends  on  L2/L1  and 
L2/T,  while  the  frequency  resolution  is  [1,  (101)  and  table  1] 

Be  -  b2  ■  O' 

(The  exact  factor  depends  on  the  particular  weighting  employed.)  The  maximum 
value  of  the  lag  weighting  length,  L2(max),  should  therefore  be  chosen  to 
achieve  the  finest  frequency  resolution  required;  smaller  values  of  L?  will 
then  result  in  coarser,  but  more  stable  estimates. 

It  can  be  seen  from  the  above  considerations  that  FFT  sizes  N1  and  N2 
in  (6)  and  (15),  respectively,  play  subordinate  roles  insofar  as  the 
fundamental  capabilities  of  this  spectral  analysis  technique  are  concerned; 
they  are  parameters  of  the  "tool"  (i.e.,  the  FFT)  being  used  to  obtain  the 
spectral  estimates,  and  must  satisfy  constraints  (11)  and  (17),  but  are 
otherwise  arbitrary.  The  fundamental  parameters  are  T,  ,  and  L2- 


SCALING 

If  we  sum  up  the  power  spectral  estimates  in  (15)  over  all  m,  we 
find  (appendix  A) 
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(19) 


m=0 


where  Pj  is  the  sample  power  of  the  total  data  record: 


T-l 


P, 


T  T 


(20) 


n=0 


Thus  if  we  adopt  the  reasonable  rule  to  maintain  the  total  sample  power,  we 


should  scale  all  the  by  the  factor 


feature 


is  incorporated  in  the  program  listings  here. 

LAG  WEIGHTINGS 

A  class  of  lag  weightings  that  encompasses  a  broad  and  useful  selection, 
including  the  rectangular,  Hanning,  Hamming,  Blackman  [3],  Harris  [19],  and 
Nuttall  [18]  windows,  is  given  by  the  form 


3 

wn  =  <*k  cos(irnk/L2)  for  jn|<  l2  . 

k=0 


(21) 


The  particular  example  we  will  concentrate  on  here  is  the"  Cl  window  in 
[18,  figure  12],  where  Cl  denotes  continuous  first  derivative,  and 


(<*k1  q  =  .355768,  .487396,  .144232,  .012604  . 


(22) 
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Since  the  lag  weighting  appears  linearly  in  the  correlation  domain,  (see 
(12)— (14)),  the  power  response  of  the  window  is  directly  proportional  to  the 
Fourier  transform  of  (21),  and  not  its  square.  Thus,  as  noted  in  [18,  (18) 
and  footnote],  the  sidelobe  level  is  half  that  depicted  in  [18,  figure  12]; 
namely,  we  have 

-46.66  dB  peak  sidelobes,  9  dB/octave  decay.  (23) 

Other  examples  are  available  in  [18]  and  can  be  used  if  deeper  sidelobes  or 
faster  decay  rates  are  required;  however,  the  main  lobe  width  must  also  be 
considered  in  these  tradeoffs. 
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CHOICE  OF  PARAMETERS 

In  this  section,  we  discuss  some  of  the  considerations  that  go  into  the 
selection  of  ,  L^,  ,  and  N^.  We  find  that  there  is  no  absolute  optimum, 

but  there  are  useful  guidelines  to  observe  in  order  to  minimize  the 
fluctuations  and  execution  time. 


STABILITY  OF  FIRST-STAGE  CORRELATION  ESTIMATES 

Because  of  the  sectioning  of  the  total  of  T  data  points  into  K  abutting 

pieces,  each  of  length  L, ,  first-stage  correlation  estimates  fa  1  in  (10) 

i  1  n 

employ 

K  L.j  =  T  data  points  for  zero  delay; 

K  (L^-l)  =  T  -  K  data  points  for  unit  delay; 

K  (L^-n)  =  T  -  n  K  data  points  for  delay  n.  (24) 

This  consideration  alone  would  say  that  to  minimize  the  loss  of  stability, 
choose  K  small,  i.e.,  segment  length  L]  large.  However,  since  lag  weighting 
(12)— (14)  will  only  use  delays  up  through  L2,  the  loss  in  stability  will  be 
relatively  insignificant  if  L2  <  L1/2,  approximately;  see  [1,  page  49, 
figure  15B]  for  specific  quantitative  results.  For  example,  the  stability  for 
window  Cl,  as  measured  in  terms  of  quality  ratio  (2)-(3),  deteriorates  by 
about  9  percent  at  L2/L^  =  1/2,  compared  to  the  absolute  optimum.  Thus  we 
will  impose  the  limitation 
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L1  >  2  l2(max)  ,  (25) 

where  lag  weighting  length  L,,(max)  corresponds  to  the  narrowest  frequency 
resolution  of  interest,  in  order  to  hold  random  fluctuations  at  an  acceptably 
low  level. 


SELECTION  OF  N., 

We  have  already  seen  in  (11)  and  figure  3  that  we  must  have  N1  >  L1  +  L2> 
in  order  to  avoid  wrap— around.  However,  there  is  no  need  for  picking  Nj 
larger,  since  the  FFTs  of  size  are  used  only  as  a  computing  shortcut  to 
getting  the  desired  correlation  of  the  averaged  L^long  data  segments,  which 
are  then  used  for  lag  weighting.  The  N^-point  forward  and  inverse  FFTs  in 
(6)  and  (10)  have  no  spectral  spacing  requirements  whatsoever,  nor  do  they 
affect  the  stability  of  the  final  spectral  estimates  {b  in  (15).  Thus  we 
will  be  interested  in  choosing  N-j  as  small  as  possible,  subject  to  both 
limitations  (11)  and  (25);  this  will  also  reduce  storage  and  increase  the 
speed  of  execution  for  each  of  these  FFTs. 


i 
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EXECUTION  TIME 


The  time  required  to  execute  one  full-precision  -point  FFT  on  the 
Hewlett-Packard  9000  Model  520  computer  has  been  determined  empirically  as 


N.|  (.14  In  N-|  +  .396)  msec  (26) 

to  a  high  degree  of  accuracy,  over  the  range  N1  =  26  =  64  to  N1  =  214  =  16384. 
The  particular  constants  in  (26)  will  change  for  a  different  computer,  but  the 
general  rule  can  be  expected  to  hold  quite  well  over  this  most  useful  range  of 
values  of  N.j .  .....  ... 

The  total  time  to  execute  all  the  K  first-stage  FFTs  in  (6)  and  the 
single  inverse  FFT  in  (10)  is  then 

( K+l )  N  (.14  In  N1  +  .396)  msec.  (27) 

To  minimize  (27),  we  should  make  and  K  small.  Since,  from  (4), 

T  =  K  L.j ,  and  since  we  want  to  use  all  available  T  data  points,  we  should 
make  large.  However,  we  are  subject  to  the  limitation  (11), 
l~|  +  (in  order  to  avoid  wrap-around),  and  we  are  interested  in 

keeping  N^  small  also,  as  noted  above,  in  addition  to  being  a  power  of  2. 
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A  reasonable  procedure  for  meeting  these  constraints  and  conflicting 
requirements  is  as  follows:  forgiven  data  length  T  and  maximum  lag  length 
L2(max)  of  interest,  choose  starter  value 

L.j  =  2  l_2(max)  (for  good  stability).  (28) 

Then  take  the  smallest  integer  n  satisfying  the  constraint 

N.|  =  2°  >  L.|  +  L2(max)  (to  avoid  wrap-around).  (29) 

Then  choose  the  largest  L1  satisfying 

L-|  <  N1  -  l_2(max)  (30) 

and  a  corresponding  K  such  that  equality 

K  L1  =  T  (to  use  all  the  data)  (31) 

is  met  (or  approximately  met). 

An  example  is  informative  at  this  point.  Suppose  T  =  1 0 i 000  and  we  want 
L2(max)  =  250.  Then  (28)-(30)  yield, in  order,  L1  =  500,  =  1024,  <  774. 

Then  if  we  take  =  774,  we  must  have  K  =  12,  yielding  K  L..  =  9288, 
which  is  significantly  below  the  allowed  value  10,000  for  this  example. 

However,  if  we  take  =  769,  then  K  =  13  yields  K  =  9997,  which  means 
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discarding  only  3  of  the  available  10,000  data  points.  Thus  a  little  juggling 
of  the  value  of  between  the  limits  given  by  (28)  and  (30),  coupled  with 
the  desire  of  meeting  (31)  as  closely  as  possible  (but  not  exceeding  T),  is 
recommended. 

This  procedure,  ( 28) — ( 31 ) ,  minimizes  execution  time  (27)  and  realizes 
near-optimum  stability  of  the  spectral  estimate.  However,  the  particular 
choices  are  not  critical.  For  example,  if  we  instead  took  =  714  and  K  =  14 
above,  we  get  K  L1  =  9996,  and  the  execution  time  (27)  increases  above  that 
for  K  =  13  by  the  factor  (14+1 )/(13+l )  =  1.071,  while  the  stability  degrades 
very  slightly;  see  [1,  figure  15B]. 

DIRECT  CALCULATION  OF  CORRELATION  * 

For  T  available  data  points,  the  number  of  multiplies  and  adds  required 
to  directly  estimate  the  correlation  at  delays  0  and  L 2  are  T  and  T  -  L2> 
respectively.  Therefore  a  total  of 


(32) 


multiplies  and  adds  are  required  for  all  delays  in  the  range  [0,  L2].  For 
the  example  above  of  T  =  10,000,  L2  =  250,  this  is  2.5E6  operations.  This 
direct  approach  takes  70  seconds  on  the  Hewlett-Packard  9000  Model  520 
computer,  versus  20  seconds  via  the  FFT  approach  presented  above.  Thus  the 
advantage  of  employing  the  FFT  technique  is  a  significant  reduction  in 
execution  time. 
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SIMULATION  RESULTS 


The  results  in  this  section  are  based  on  the  two  program  listings  in 
appendix  B  for  complex  and  real  data,  respectively.  An  explanation  of  their 
use  is  given  at  the  beginning  of  that  appendix. 


SINGLE  TONE,  NOISE  FREE 


The  first  example  is  a  noise-free  pure  complex  tone  at  zero  frequency. 


with  A  =  1  Qprnnrl  T  =  10  000  data  nnintc  I  =  7AQ  K  =  17  niarpc 
"  •  ’  - *  *  '  ~  » - - r  ~  *  ‘--j  .%  ,w  r  ■  -  -  -  ^  , 


N^  =  1024,  and  the  Cl  weighting  given  in  (21 )— (23) .  The  resulting  spectral 
estimate  for  L^  =  250,  N^  =  16384,  obtained  via  the  programs*  in  appendix 
B,  is  given  in  figures  4  and  5,  for  frequencies  f  =  m/N2  in  the  range 
[0,  8/L2];  the  first  null  is  at  frequency  4/(2L2)  =  2/L2.  (Compare  this 
spectrum  with  the  ideal  result  in  [18,  figure  12].)  The  difference  in  figures 
4  and  5  is  that  the  dB  measure 


1°  log  |Bj 


(33) 


*Observe  that  constraints  (11),  (17),  and  ( 28)— ( 31 )  have  been  observed,  both 
in  this  example  and  in  the  programs. 
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Figure  4.  Pure  Tone  Analysis  for  N2=  16384; 
Negative  Lobes  Also 


Figure  5.  Pure  Tone  Analysis  for  N2 =  16384; 
Positive  Lobes  Only 
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Is  plotted  on  the  ordinate  in  figure  4,  whereas  the  negative  lobes  (Bm<0) 

are  not  plotted  in  figure  5.  The  latter  procedure,  not  (33),  is  the  correct 

one  and  is  adopted  henceforth.  The  reason  for  the  negative  lobes  in  the 

spectral  estimate  is  that  the  lag  weighting  and  lag  window  appear  linearly,  not 

quadratically,  in  the  equations;  see  (22)  et  seq.  and  [1,  ( 50) -( 51 ) ,  (97 ) -( 1 00) ] . 

The  presence  of  negative  spectral  estimates  for  some  frequencies,  i.e.,  B  <0 

m 

for  some  m,  can  actually  be  useful  as  an  indicator  of  the  presence  of  strong 
narrowband  components,  as  noted  by  Blackman  and  Tukey  [3,  pages  13,  92,  115]. 

Of  course,  the  addition  of  noise  would  fill  in  the  deep  valleys  in  figures  4 
and  5. 


The  danger  of  using  too  small  a  value  of  for  the  second -s  tage  FFT  is 
depicted  in  figures  6  and  7  for  N2  =  1024  and  512,  respectively.  The  only 
difference  in  figures  5,  6,  and  7  is  the  value  of  ^(generally,  any 
parameters  not  mentioned  are  maintained  at  the  same  values  as  for  the  previous 
figure).  The  frequency  spacing,  (16),  becomes  progressively  coarser,  to  the 
point  that  the  linear  Interpolation  between  frequency  components  (15), 
employed  in  all  the  figures,  can  lead  to  some  questionable  conclusions;  hence 
we  should  observe  requirement  (17)  on  N^.  Since  only  one  FFT  of  size  N2 
need  be  performed,  for  each  choice  of  lag  length  l_2,  this  is  not  a  severe 
computational  load;  the  bulk  of  the  first-stage  processing,  using  several 
small-size  (N^ )  FFTs  is  done  only  once  and  saved  for  as  many  second-stage 
spectral  estimates  of  different  resolutions  as  desired.  Of  course,  for  the 
coarser  frequency  resolutions,  L2  is  smaller,  thereby  alleviating  the 
requirement  (17)  on  N  . 
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Figure  7.  Pure  Tone  Analysis  for  N2 =  512 
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Since  there  is  no  additive  noise  in  this  particular  example,  there  is 
really  no  need  to  observe  constraint  (28)  here.  If  we  keep  l_2  fixed  at  250, 
we  can  decrease  L1 ,  still  subject  to  (12) -(14) ,  and  also  decrease  N1 , 
subject  to  (11),  and  realize  the  same  spectral  estimate  as  in  figure  5.  An 
example  for  L-j  =  256,  K  =  39,  K  L1  =  9984,  N1  =  512,  N2  =  16384,  yielded 
a  result  indistinguishable  from  figure  5,  and  is  not  plotted  here.  This 
confirms  the  earlier  conclusion  that  the  effective  window  depends  only  on  the 

desirable  weighting  jv/ in  (13)— (14) ,  and  not  on  the  particular  choices  of 
L,.  K,  Nr 

In  figure  8,  the  frequency  of  the  pure  tone  is  changed  from  zero  to  1/4 
Hz,  and  the  entire  spectral  estimate  over  frequency  range  (-  1/2,  1/2)  Hz  is 
plotted,  with  N2  returned  to  value  16384.  The  rapid  decay  and  deep 
sidelobes,  (23),  associated  with  the  Cl  weighting  in  (21 ) — ( 22 )  are  quite 
evident.  The  darkened  portion  of  the  plot  is  caused  by  the  detailed  sidelobe 
structure  of  the  Cl  window. 

On  the  other  hand,  if  we  simply  eliminate  the  division  of  the  first-stage 
correlation  estimates  by  the  autocorrelation  of  the  temporal  weighting,  i.e., 
eliminate  the  triangular  1  -  n/L1  term  in  (13)  and  (14),  then  the  spectral 
estimate  in  figure  9  is  obtained.  Despite  the  retention  of  the  desirable 
weighting  in  (13)  and  (14),  there  is  a  significant  fill-in  of  the  deep  valleys 
and  a  less  rapid  decay  of  the  estimate  in  figure  9. 
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Frequency  f 

Figure  8.  Pure  Tone  at  Frequency  .25  Hz; 

Lag-Reshap i ng 


-.5  -.25  0  .25  .5 

Frequency  f 


Figure  9.  Pure  Tone  at  Frequency  .25  Hz; 

No  Lag-Reshaping 
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TWO  SEPARATED  TONES  IN  WHITE  NOISE 


A  similar  comparison  of  the  effect  due  to  lack  of  lag-reshaping,  i.e., 
removal  of  division  by  1  -  n/L^ ,  is  displayed  in  figures  10  and  11  for  two 
wel 1 -separated  pure  tones  plus  white  noise.  The  strong  tone  at  frequency  .25 
Hz  has  a  power  level  47.8  dB  stronger  than  the  total  white  noise  power,  while 
the  weak  tone  at  frequency  .0625  Hz  (indicated  by  an  arrow)  has  a  power  level 
-12.2  dB  relative  to  the  noise.  Thus  the  ratio  of  signal  powers  is  60  dB. 
Whereas  there  is  an  indication  of  the  weak  tone  in  figure  10,  at  the  correct 
level,  there  is  none  in  figure  11,  because  of  the  poor  sidelobe  structure  in 
the  latter  procedure.  In  fact,  the  sidelobes  dominate  the  noise  spectrum 
completely  in  figure  11. 

In  figures  12  and  13,  the  exact  procedures  are  repeated  for  the  same 
parameter  values,  except  that  the  data  are  restricted  to  be  real  and  the 
spectrum  is  only  plotted  over  frequency  range  (0,  1/2)  Hz.  (The  second 
program  listing  in  appendix  B,  for  real  data,  was  employed  for  these  two 
figures.)  Again,  the  weak  signal  is  indicated  only  for  the  case  of 
lag-reshaping  in  figure  12. 
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Frequency  f 

Figure  10.  Two  Complex  Tones  plus  White  Noise; 

Lag-Reshap i ng 


Figure  11.  Two  Complex  Tones  plus  White  Noise; 

No  Lag-Reshaping 
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Figure  12.  Two  Real  Tones  plus  White  Noise; 

Lag-Reshap i ng 


Figure  13.  Two  Real  Tones  plus  White  Noise; 

No  Lag-Reshaping 


30 


TR  7511 


TWO  CLOSE  TONES  IN  WHITE  NOISE 


In  figures  14  and  15,  a  strong  complex  tone  at  frequency  .25  Hz  has  a 
power  level  17.8  dB  stronger  than  the  total  white  noise  power,  while  the  weak 
tone  at  frequency  .265  Hz  has  a  power  level  of  -12.2  dB  relative  to  the 
noise.  Thus  the  ratio  of  signal  powers  is  30  dB.  Whereas  the  lag-reshaping 
spectral  estimate  in  figure  14  succeeds  in  resolving  the  two  close  tones  at 
frequency  separation  .015  Hz,  the  result  in  figure  15,  for  no  lag-reshaping, 
indicates  only  one  of  the  tonals.  It  may  be  observed  that  the  detailed 
wiggles  in  the  two  spectral  estimates  are  virtually  identical,  except  for 
frequencies  near  the  tone  locations. 

Finally,  in  figures  16  and  17,  the  effect  of  decreasing  lag  weighting 
length  L2  is  demonstrated.  Specifically,  figures  14,  16,  17  all  employ 
lag-reshaping,  the  only  difference  in  the  three  figures  being  that  L?  =  250, 
200,  150,  respectively.  Whereas  the  weak  close  tonal  Is  wel 1 -resol ved  in 
figure  14,  it  is  lost  in  figure  17;  the  estimate  in  figure  16  is  Intermediate 
and  on  the  border  of  being  resolved.  The  familiar  tradeoff  of  resolution 
versus  stability  is  well  demonstrated  in  figures  14,  16,  17.  Namely,  as  the 
resolution  degrades  (mainlobe  widens),  the  fluctuations  in  the  spectral 
estimate  decrease;  figures  16  and  17  are  progressively  smoother  versions  of 
figure  14.  Which  case  to  prefer  depends  on  the  particular  situation  under 
investigation,  including  factors  such  as  the  proximity  of  tones  with  widely 
different  strengths,  and  on  the  particular  colored  noise  spectrum 
encountered.  Quantitative  evaluation  of  the  stability  Is  considered  in  the 
next  section. 
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-.5  -.25  0  .25  .5 


Frenuencv  f 

- ,  -  ~  y 

Figure  14.  Two  Close  Tones  plus  White  Noise; 
Lag-Reshaping,  L2 =  250 


-.5  -.25  0  .25  .5 

Frequency  f 

Figure  15.  Two  Close  Tones  plus  White  Noise; 
No  Lag-Reshaping,  Lz =  250 
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Figure  16.  Two  Close  Tones  plus  White  Noise; 
Lag-Reshaping,  L2 =  200 


-.5  -.25  0  .25  .5 

Frequency  f 

Figure  17.  Two  Close  Tones  plus  White  Noise; 
Lag-Reshaping,  L2 =  150 
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VARIANCE  OF  SPECTRAL  ESTIMATE 


THEORETICAL  RESULTS 

The  variance  of  the  spectral  estimate  obtained  via  the  generalized 
technique  Including  temporal  and  lag  weighting  was  derived  and  evaluated  in 
[1,  pages  41-57],  In  particular,  for  abutting  rectangular  temporal  weighting, 
the  NQR,  defined  in  (2)  et  seq.  above,  is  given  by  [1,  (138)]  as 


NQR  = 


(34) 


0 


in  the  case  of  lag-reshaping,  where  w^Cr)  is  the  continuous  analog  of  the 


(d) 

discrete  desirable  weighting  wn  employed  here  in  (13)  and  (14). 


For  the  class  of  continuous  weightings  (compare  (21)) 


(35) 


the  NQR  in  (34)  becomes 


Z 


0  -  x  Lg/L-j ) 


-1 


NQR  = 


(36) 
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which  depends  only  on  the  ratio  L^/L-j  and  coefficients  This 

quantity  was  plotted  in  [1,  page  49,  figure  15B]  for  various  windows.  The 
special  case  of  the  Cl  window,  (22),  is  replotted  here  in  figures  18  and  19  as 
the  solid  curve. 

SIMULATION  RESULTS 


In  order  to  corroborate  these  theoretical  results,  white  Gaussian  noise 
was  generated  for  25  independent  trials,  each  of  length  T  =  10000  samples. 
The  data  were  spectral  analyzed  for  lag  length  L^  =  1000  and  for  a  variety 
of  temporal  lengths  L^ .  The  results  for  the  NQR,  determined  by  using  alj_ 
frequency  bins  on  all  the  trials,  are  shown  as  crosses  in  figure  18.  The 
optimum  quality  ratio,  as  given  by  (3),  is 


2  *  1000  *  .2558 

10000 


.0512  , 


(37) 


where  we  employed  [1,  (11),  (101),  and  Table  1].  Keeping  L2  fixed  is 
tantamount  to  holding  the  frequency  resolution  constant,  while  varying  L^ 
corresponds  to  different  segment  lengths  (all  subject  to  L1  >  L2,  of 
course) . 


Results  for  L2  =  250  are  plotted  in  figure  19,  based  now  on  at  least 
100  trials  for  each  value  of  L1  considered;  now  the  optimum  Q  =  .0128.  The 
simulation  results  for  white  Gaussian  noise  slightly  overestimate  the 
theoretical  results  in  figure  19,  whereas  the  converse  is  true  in  figure  18. 
Also  added  in  figure  19  are  simulation  results  for  white  input  data  which  are 
not  Gaussian,  but  rather  have  a  flat  probability  density  function  over  their 
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nonzero  extent;  see  the  circles,  labeled  "rectangular  data,"  in  figure  19. 

The  discrepancy  between  the  theoretical  results  for  Gaussian  data  and  the 
simulation  results  for  rectangular  data  is  believed  to  be  real  (although  it  is 
only  approximately  a  2  percent  difference),  because  these  results  were  based 
on  500  trials  for  each  value  of  Lj .  In  any  event,  the  theoretical  results 
predict  the  stability  very  well,  over  the  full  range  of  L^/L^ .  At 
l^/L-!  =  .5,  for  example,  the  loss  in  stability  is  only  about  8  to  9 
percent,  as  measured  by  the  quality  ratio  (2). 


APPROXIMATION  TO  NQR 

For  small  l^/L-j ,  the  NQR  in  (34)  behaves  approximately  as 

„n  'o  +  I!  (W  <W2  +  I3  <VLi>3 

NIJK  =  - 7 -  , 

0 

where  we  define  the  dimensionless  constants 
1 

In  =  f  dx  wjj(L2x)  xn  for  n  =  0,  1 ,  2,  3  . 

0 

For  the  class  of  continuous  weightings  specified  in  (35),  there  follows 


(38) 


(39) 


2  12  2  2 
*0  =  a0  +  2  *al  *  a2  +  a3*  ’ 


I  =XI  -  —  T 

h  2  A0  21’ 

ir 


T  _  1  T  4.  -1-  (  2  4.  1  3  1  2,4_t 

l2  3  !0  .  2  (al  4  a2  9  a3^  2  T2  ’ 

4ir  v 


(40) 
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where 

T1  =  a0  al  +  9  a0  a3  +  9  al  a2  +  25  a2  a3  ’ 

T2  =  T1  “  4  a0  a2  '  32  al  a3  '  (41 > 

Equations  ( 38) -(41 )  enable  a  simple  approximate  calculation  of  the  NQR  for  any 
weighting  in  class  (35). 


When  applied 

to  the  Cl , 

C3,  C5 

windows  in  [18, 

figures  10,  11 ,  12],  the 

numerical  results 

for  (Vo! 

[3  are, 
'l 

respectively, 

.16114, 

.04036, 

.01273 

for  window 

Cl; 

.15362, 

.03664, 

.01100 

for  window 

C3; 

.14168, 

.03111 , 

.00859 

for  window 

C5.  (42) 

Substitution  in  (38)  yields  results  which  agree  very  well  with  figures  18  and 
19  and,  more  generally,  with  [1,  figure  158]. 
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SUMMARY 


Programs  for  spectral  estimation  of  complex  data  as  well  as  real  data 

have  been  given  and  exercised  under  a  wide  range  of  parameter  choices.  They 

allow  the  data  to  be  generated  or  made  available  in  abutting  blocks  of  L1 

data  points  at  a  time,  and  thay  employ  no  temporal  weighting  at  all.  The  user 

must  replace  subroutine  SUB  Data  by  his  own  data  input  routine.  The 

subroutine  SUB  Fftl4z  listed  in  appendix  B  employs  zero-subscripted  arrays,  as 

14 

encountered  directly  in  the  equations,  and  can  handle  sizes  up  to  2 
16384. 


The  benefits  to  be  accrued  from  lag-reshaping  have  been  demonstrated  by 
simulation,  for  a  variety  of  situations  including  multiple  tones  in  noise. 

The  ability  to  obtain  several  spectral  estimates  with  different  frequency 
resolutions  has  been  Incorporated  in  the  programs  in  such  a  way  that  the  bulk 
of  the  first-stage  spectral  estimates  do  not  have  to  be  recalculated.  This 
enables  the  user  to  make  his  own  decisions  on  resolution  versus  stability, 
without  having  to  do  extensive  computations  repeatedly. 
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APPENDIX  A.  SUM  OF  POWER  SPECTRAL  ESTIMATES 


The  power  spectral  estimate  Bm  is  given  by  (15)  for  0  <  m  <  N?  -  1. 
The  sum  over  all  m  is 


N  -1 


B  =  N,  b„ 
m  2o 


N~  w<d>  a 
2  o  o 


m=0 


N  w(d) 
n2  o 


Nrl 


N,  -1  K 


A_  = 


N  w(d) 
n2  wo 


m=0 


,00 

'm 


m=0  k=l 


K  N  -1 


K  L  -1 


N  w(d) 
2  0 


2^  Rk>r  -  «id>  i  i 


k=l  m=0 


k=l  n=0 


,00 


T-l 


-  M1  N2  “o'”  2  M2  • 


n=0 


(A-l ) 


by  use  of  (14),  (10),  (9),  (6),  (5),  (4),  in  order.  In  terms  of  the  sample 
power  of  the  total  data  record, 


T-l 


PT  ■  T  2  hi2  • 

n=0 


(A-2) 


the  sum  in  (A-l)  becomes 


N  -1 


m=0 


8  =  N.  N,  K  L.  wid)  PT 

m  12  1  o 


(A-3) 
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APPENDIX  B.  PROGRAMS  FOR  LAG-RESHAPING  SPECTRAL  ANALYSIS 


There  are  two  main  programs  in  this  appendix,  written  for  the 
Hewlett-Packard  9000  Model  520  computer.  The  first  is  for  complex  data,  while 
the  latter  accommodates  only  real  data.  The  following  explanation  is  keyed  to 
the  complex  data  program,  but  applies  directly  to  the  real  data  program  as 
well. 


Input  parameters  T,  L1 ,  L2(max),  N1  are  required  at  the  top  of  the 
program  in  lines  20-50.  The  lag  weighting  coefficients  of  interest  are 
entered  in  line  60.  Constraint  (29)  is  enforced  in  lines  100-130.  The  number 
of  pieces,  K,  and  the  number  of  data  points  used,  K  L1 ,  are  computed  in 
lines  160  and  180.  The  input  data  are  entered  via  SUB  Data  in  line  280,  L1 

data  points  at  a  time;  that  is,  data  values  jx^^J  for  0  <  n  <  L1  -  1 ,  as 
given  by  (5),  are  accessed  by  the  CALL  in  line  280  with  arguments  Ks,  LI  (= 
piece  k,  segment  length  L^ ) .  Division  by  the  autocorrelation  of  the 
temporal  weighting  is  accomplished  once,  in  loop  380-420. 

At  this  point,  the  first-stage  calculations  are  complete;  the  correlation 
domain  quantities  are  stored  in  arrays  Xa,  Ya.  Input  parameters  L2>  N2 
must  now  be  entered  in  lines  440,  450.  Constraints  (14)  and  (11)  are  enforced 
in  lines  460-510.  The  lag  weighting  is  applied  in  loop  620-680,  taking  figure 
3  into  account,  while  the  scaling  considerations  of  ( 1 9 ) -( 20)  imply  line  730. 
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Finally,  the  spectral  estimate  is  plotted  in  dB  over  the  frequency  range 
(-1/2,  1/2)  Hz  in  lines  770-860.  (For  real  data,  the  spectrum  is  only  plotted 
over  the  range  (0,  1/2)  Hz,  since  it  is  symmetric  about  the  origin.) 


The  FFT  subroutine  in  lines  900  -  1710  uses  zero-subscripted  arrays,  just 

as  encountered  in  the  usual  mathematical  definition;  it  can  handle  FFT  sizes 

14 

up  to  2  =  16384.  The  subroutine  SUB  Data  in  lines  1730  -  1900  must  be 

replaced  by  the  user;  however,  notice  should  be  taken  of  lines  1810  -  1830  and 
the  storage  locations  in  lines  1860  -  1870,  in  order  to  ensure  that  data 
points  at  a  time  are  properly  stored  in  locations  0  through  -  1 . 


10 
c-  y 
3Q 
40 
5G 
60 
70 
80 
90 
1  0  0 
110 
120 
1  30 
'  140 
150 
1  60 
170 
1  8  0 
190 
200 


LAG  RE-SHAPING 
1=18080 
LI =769 
L 2m ax =2 50 
N 1=1 024 


SPECTRAL  ANALYSIS  FOR  COMPLEX  DATA 

!  TOTAL  NUMBER  OF  BATA  POINTS 
!  NUMBER  OF  BATA  POINTS  PER  PIECE 
!  MAXIMUM  LAG  WEIGHTING  LENGTH 
!  SIZE  OF  FIRST-STAGE  FFT  N 1 >=L 1 +L2max 


BATA  •.  355768,  .487396,  ,  144232 ,  .012604  !  LAG . WEIGHT ING  COEFFICIENT 

T  R  6239,  p  ag  e  1  6 ,  figure  12:  -46.66  d  B  *s  i  del  o  h  e  s ,  9  d  E  /  o  c  t  av  e  d  e  c  ay 
D  I M  X  <0:1  6383  >  ,  Y (0M 6383 >  ,  Xa<  0 :  8 1 9 1  )  ,  Ya< 0 :  8 191  > 

DOUBLE  T,  LI  ,  L2max,  N1  ,  Nlrn,  K,  Ks,  Ns,  L2,  N2, N2m  !  INTEGERS 

IF  N 1 >  =  L  1  +  L  2  rri  ax  THEN  140' 

PRINT  "HI  INCREASED,  SO  AS  TO  EE  AT  LEAST  Ll+L2max" 

N  1  m= I  NT  < LOG < L 1 +L2max  > /LOG <  2) >  + 1 


N 1  =  2  •A-  N  1  in 


N 1  m  =  N 1  - 1 

R  E  D  I M  X  <  0  :  N  1  m  >  ,  Y  <  0 :  N 1  m  )  ,  X  a  <  0  :  N  1  m  >  ,  Y  a  <  0  :  N  1  m ') 

K= I  NT  <  T/L 1 >  !  NUMBER  OF  DATA  SEGMENTS 

PRINT  "T  =" ;T;"  LI  =  "  ;  L 1 ;  "  L2max  = "  ;  L2max  ;  M  N1  =  " ;  N 1 5  11  K  =!l 

PRINT  "USED" 5  K*L 1 ; “OF  THE  TOTAL  NUMBER  OF  DATA  POINTS, " ; T 

READ  A 0 , A  1 , A 2 , A3 

A  = 1 . / < A 0 + A 1 +  A 2 + A 3 )  !  NORMALIZE  LAG  WEIGHTS 


210  A  0  =  A  0  *  A 


220  A  1 = A  1  * A 


230  A2=A2*A 

240  A  3 = A  3  *  A 

250  MAT  Xa=<0. > 

260  MAT  Y a= < 0 .  > 

270  FOR  K s  = 1  TO  K 

280  CALL  D at  a  < K  s , L 1 , X  <  *  >  ,  Y <  *  >  >  !  LI  DATA  POINTS  AT  A  TIME 

290  FOR  Ns=L 1  TO  Him 
3O0  X  <  Ns  >  =  Y  <  Ns  >  =  0 . 

310  NEXT  Ns 
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328 

CALL  Fftl4z<Nl,X<*>, 

Y  (  * >  ) 

!  INTO  FREQUENCY  DOMAIN 

330 

FOR  N s - O  TO  Him 

340 

X  a  < Ns > - X a  < Ns > +X ( Ns  > * 

X  <  Ns  >  +  Y 

< Ns > *Y C Ns >  !  POWER  AVERAGING 

350 

NEXT  Ns 

360 

NEXT  Ks 

370 

CALL  Fft 14z<Nl,Xa<*> 

,  Y  a  *  >  > 

!  INTO  LAG  DOMAIN 

330 

FOR  N s  =  0  TO  Ll-1 

390 

A  =  L 1 / <  L 1 -Ns> 

400 

X  a  <  N  s  )  =  X  a  <  N  s  >  *  A 

!  DIVIDE  BY  CORRELATION 

410 

Y  a  <  N  s  >  =  -  Y  a  <  N  s  )  #  A 

!  OF  TEMPORAL  WEIGHTING 

420 

NEXT  Ns 

430 

BEEP 

440 

INPUT  " L2  =" ,L2 

!  LAG  WEIGHTING  LENGTH  L2CL1,  L2<= 

458 

INPUT  11 N  2  =  ",N2 

!  SI 

ZE  OF  SECOND-STAGE  FFT  N2  OL2 

460 

IF  L2CL1  THEN  490 

470 

PRINT  " L2  REDUCED .TO 

Ll-1, 

TO  AVOID  DIVISION  BY  ZERO" 

480  • 

L2-L 1 - 1 

490 

IF  L2<  =L2max  THEN  52 

0 

500 

PRINT  " L2  REDUCED  TO 

L  2  m  ax  , 

TO  AVOID  WRAP-AROUND  REGION" 

510 

L  2  =  L  2  m  ax 

520 

IF  N  2  >  L  2  +  L  2  THEN  560 

530 

PRINT  " N 2  INCREASED, 

SO  AS 

TO  BE  GREATER  THAN  L2*2" 

540 

N2m= I NT  <  LOG<  L2+L2+ 1 > 

/  L  0  G  <  2  > 

>  +  l 

550 

N  2 = 2  N  2  m 

560 

PRINT  "L2  =st ;  L2  ;  " 

N2  =" 

5  N2 

570 

N  2  m  =  N  2  - 1 

580 

REDIM  X<0:  N2rn>  ,  Y<0: N2m> 

590 

A  =  P  I  /  L  2 

600 

X  C  0  >  =  X  a  (  0  ) 

610 

Y  <  0  >  =  0 . 

620  FOR  Ns  = 1  TO  L2 

630  B  =  N  s  *  A 

640  W=A0+A1*COS<B)+A2*COS< B  +  B > +R3*C0S < B  +  B  +  B > 

650  X <  Ns >  =X < N2-Ns >=Xa< Ns > * W  !  LAG  WE  I GHT I NG 

660  Y  ( Ns  >  =B=*Ya<  Ns  >*W 

670  Y<N2-Ns)=-B 

680  NEXT  Ns 

690  FOR  N s  =  L 2  +  1  TO  N2rr.-L2 

700  X  <  Ns ) -Y < Ns  >  =  0 . 

710  NEXT  Ns 

720  CALL  Fftl4z<N2,X<*>,  Y<*)>  !  INTO  FREQUENCY  DOMAIN 

730  MAT  X  =  X s < FLT < K ) *L 1 *N 1 *N2 > 

740  Bi g=MAX<X<*> > 

750  D  b = L  G  T <  B  i g  )  *  1 0 . 

760  PRINT  "Big  =“;Big; "  dB  max  =" ;  Db 

770  PLOTTER  IS  "GRAPHICS" 

73 0  GRAPHICS  ON 

790  WINDOW  -N2X2,  H2/‘2, -70.  ,  0. 

800  GRID  N 2 /■  3 ,  10. 

310  FOR  Ns=-N2/2  Tu  N2---2 

320  Ks=Ns  MODULO  N2 

830  W=MAX<X<Ks> , 1 . E-20> 

840  PLOT  Ns , LGT < W ) *  1 0 . -Db 

350  NEXT  Ns 

860  PENUP 

870  GOTO  440 

330  END 

890  ! 
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900  SUE  Fft  14z<  DOUBLE  N, REAL  X(0,YC*))  !  .  N  =  2A  I NTEGERC  =  2* 1  4=  1  6384 ;  Q  SUBS 

910  DOUBLE  Ml  ,  N2,  N3,  N4 ,  Log2n,  J,  K  !  INTEGERS  <  2--31  =  2,147,483,643 

928  DOUBLE  II, 12, 13, 14,15',  16,17,18,19,110,111,112,113,114 
930  N  1  -  H  ■■■'  4 

940  ALLOCATE  C < 0 : H 1 >, DOUBLE  L<0: 13)  !  QUARTER-COSINE  TABLE  IN  C<*> 

950  H=(FI+FT)/N 

960  FOR  J=0  TO  Ml 

970  C< J)=CGS< A*J) 

980  NEXT  J 

990  N2=N1+1 

1000  N  3  =  N  2  + 1 

1010  N  4  =  N 3  +  N 1 

1020  L  o  g  2  n  =  1 . 4 4 2 7 L 0 G  ( N  > 

1030  FOR  11  =  1  TO  L o g 2 n 

1040  I  2  =  2  A  C  L  o  g  2  n  -  I  1 ) 

1050  13=12+12’ 

1060  I  4  =  N I  3 

1070  FOR  15=1  TO  12 

1080  1 6=  C I  5- D  *  1 4+ 1 

1090  IF  I  6 <  =  N 2  THEN  1130 

1100  A  1  =  -  C  C  N  4  - 1  6  -  1 > 

1110  A2=-C< I6-N1-1  > 

1120  GOTO  1150 

1 130  A 1 =  C  < 16-1  > 

114S  A 2  =  - C  <  N 3 - I  6 -  1  > 

1150  FOR  17  =  0  TO  N - 1 3  STEP  13 

1160  13=17+15-1 

1170  19=13+12 

1180  A 3  =  X  < I  3  > - X ( 19) 

1190  A 4  =  V  < 1 8 ) - Y  < I  9  > 

1200  XC I8)=X( I8)+X(I9)  v 

1210  Y(I8)=Y(I8)+Y(I9) 


1220 

X<I9>=A1*A3-I 

A2*A4 

1230 

Y<I9>=A1*A4+A2*A3 

1240 

NEXT 

17 

1250 

NEXT 

15 

1260 

NEXT 

I  1 

1270 

1 1  =  L  o  g  2  n  +  1 

1280 

FOR 

12=1  TO 

14 

1290 

LU2 

-1  )  =  1 

1  300 

IF  I 2  >  L  o  g  2  n 

THEN  1320 

1310 

L(  12 

-1)=2M  11-12) 

1320 

NEXT 

12 

1330 

K  =  0 

1340 

FOR 

11=1  TO 

L  <1 3  > 

1350 

FOR 

12=11  TO 

LC12)  STEP  L(13) 

1360 

FOR 

13=12  TO 

L ( 1 1 )  STEP  L(12) 

1370 

FOR 

14=13  TO 

L < 1 0 )  STEP  L (  1  1  ) 

1380 

FOR 

15=14  TO 

LC 9 )  STEP  L < 1 0 ) 

1390 

FOR 

16=15  TO 

L < 8 >  STEP  L < 9 ) 

1400 

FOR 

17=16  TO 

L< 7 )  STEP  L < 3 ) 

1410 

FOR 

18=17  TO 

L < 6 >  STEP  LC7) 

1  420 

FOR 

19=13  TO 

L<5)  STEP  L  < 6 > 

1  430 

FOR 

110=19  TO  L<4)  STEP  LC5) 

1440 

FOR 

111=110 

TO  Lc.3)  STEP  LC4) 

1450 

FOR 

112=111 

TO  L  < 2 )  STEP  L  <  3  > 

1460 

FOR 

113=112 

TO  LC1)  STEP  L C 2 ) 

1470 

FOR 

114=113 

TO  L < 0 >  STEP  LCD 
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1480 

1490 

1500 

1510 

1520 

1530 

1540 

1550 

1560 

1570 

1580 

1590 

1600 

1610 

1620 

1630 

1640 

1650 

1  660 

1670 

1630 

1690 

1700 

1710 

1720 

1730 

1740 

1750 

1760 

1770 

1780 

1790 

1300 

1310 

1320 

1830 

1340 

1350 

I860 

1870 

1330 

1390 

1900 
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J  = I  14-1 

IF  fOJ  THEN  1560 
A  =  X  GO 
XGO=X<  J) 

X<  JV=A 
fl  =  YGO 
YGO=Y<  J> 

Y  (  J  )  =  fi 
K  =  K+ 1 
NEXT  114 
NEXT  113 
NEXT  112 
NEXT  Ill 
NEXT  110 
NEXT  19 
NEXT  13 
NEXT  17 
NEXT  16 
NEXT  15 
NEXT  14 
NEXT  13 
NEXT  12 
NEXT  II 
SUBEND 


SUB  Bat a< DOUBLE  Ks, LI , RERL  X<*>,Y<*>> 
DOUBLE  J  1  ,  J  2  , J 


W 1 =2 . *PI*. 25 
W  2  =  2  . *  P  I  *  . 265 
T  1  =  .  3  1 
T2=. 77 
Db2=-30 . 


TONE  FREQUENCY  1 
TONE  FREQUENCY  2 
TONE  PHASE  1 
TONE  PHASE  2 
RELATIVE  TONE  STRENGTH 


A  2  = 1 0 , A<Db2/20. > 

J  2  =  K  s  *  L  1  -  1 
J 1 = J2-L 1+1 
FOR  J=J1  TO  J2 
P1=W1*J+T1 
P  2  =  W  2  *  J  +  T  2 

XC J-Jl )=C0S<P1 >+fl2*C0S<P2)+CRND-. 5>*. 316 
Y  < J- J 1 > =S I N  < P 1 >  +  A2*S I N  C P2 )  +  <  RND- . 5  >  * . 3 1 6 
!  CRND-. 5>  +  l < RND- . 5 >  has  power  1/12+1/12  =  -7.73  dB 
NEXT  J 
SUBEND 
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10 
20 
30 
40 
50 
60 
70 
80 
90 
100 
1  1  0 
120 
130 
140 
150 
160 
170 
1 8  0 
190 
200 
210 
220 
230 
240 
258 
260 
270 
230 
290 
300 
310 
320 
330 
340 
350 
360 
370 
380 
390 
400 
410 
420 
4  3  0 
440 
450 
460 
470 
480 
490 
500 
510 
520 
530 
540 
550 
560 
570 


LAG  RE-SHAPING  SPECTRAL 
T  =  1  0  0  0  0  ! 

Li  =769  ! 

L2rnax=250  ! 

N 1  =  1 0  2  4  ! 

DATA  . 355768 , . 487396, . 144232, 
TR  62 3 9 ,  p ag e  16,  f i g u r e  12: 


ANALYSIS  FOR  REAL  DATA 
TOTAL  NUMBER  OF  DATA  POINTS 
NUMBER  OF  DATA  POINTS  PER'  PIECE 
MAXIMUM  LAG  WEIGHTING  LENGTH 
SIZE  OF  FIRST-STAGE  FFT  N  1  >  =  L 1  +  L  2 m ax 

012604  !  LAG  WEIGHTING  COEFFICIENTS 

46.66  d  E  side!  o  b  ■=■  s ,  9  d  B  /  o  c  t  av  e  d  e  c  ay 


DIM  X < 0 :  163 8 3 > ■ 


:  0 : 1 6383 > , X a < 8:8191) 


i 


DOUBLE  T  ,  L  1  ,  L2max ,  N  1 ,  N  1  m ,  K ,  Ks ,  Ns  ,  L2  ,  N2 ,  N2rn 
IF  HI >  =  L 1 +  L 2 m ax  THEN  148 

PRINT  “Nl  INCREASED,  SO  AS  TO  EE  AT  LEAST  Ll+L2max" 

N 1 rn  = I  NT  < LOG <  L 1 +L2max > /LOG < 2 ) >  + 1 

Nl=2-Nlm 


INTEGERS 


N 1  rn = N 1  - 1 

RED  I M  X  <  8 :  N  1  m  ) ,  Y  <  8 :  N 1  rn  >  ,  Xa<  0  :  N  1  rn  ) 

K= I  NT  <  T/L 1 >  .  !  NUMBER  OF  DATA  SEGMENTS 

PRINT  "T  =  "  ;T;M  LI  = 11  ;  L 1  ;  "  L2rnax  =  n  ;  L2rr.ax ;  “  Nl  =  11 ;  N  1 ;  11  K  =  ";* 

PRINT  "USED" ; K*L1 5 "OF  THE  TOTAL  NUMBER  OF  DATA  POINTS, " ; T 

READ  A 0 , A  1 , A 2 , A3 

A=1 . /( A0+A1+ A2+A3)  !  NORMALIZE  LAG  WEIGHTS 

A  0  =  A  8  *  A 

A 1  =  A  1  *  A 

A2= A2* A 

A  3  =  A  3  *  A 

MAT  Xa-<8. ) 

FOR  Ks= 1  TO  K 

CALL  D  at  a  <  K s  ,  L 1 , X < * ) )  !  LI  DATA  POINTS  AT  A  TIME 

FOR  N  s  =  L  1  TO  N 1  rn 
X<  Ns ) =0 . 

NEXT  Ns 
MAT  Y =  <  8 .  ) 

CALL  F f t 1 4 z < N 1 , X < * ) , Y < * > >  !  INTO  FREQUENCY  DOMAIN 

FOR  N s  =  8  TO  N 1  rn 


X a < Ns ) =  X a < Ns) +X < Ns )*X < Ns ) + Y < Ns ) #Y < Ns >  !  POWER  AVER AG  I NG 

NEXT  Ns 
NEXT  Ks 
MAT  Y=<0. ) 

CALL  Ff t 14z<Nl  ,  Xa<*)  ,  Y<*> )  !  INTO  LAG  DOMAIN 

FOR  N s  =  0  TO  L 1 - 1 

X a < Ns ) =Xa< Ns ) *L 1 / < L 1 -Ns  >  !  D I V I DE  BY  CORRELAT I  ON 

NEXT  Ns  !  OF  TEMPORAL  WEIGHTING 

BEEP 

INPUT  “ L2  ="  ,L2  !  LAG  WEIGHTING  LENGTH  L2CL1,  L20N1-L1 

INPUT  " N2  = “ , N 2  !  SIZE  OF  SECOND-STAGE  FFT  N2  <  >L2*10  IS  GOOD 

IF  L 2 < L 1  THEN  480 

PRINT  “L2  REDUCED  TO  Ll-1 ,  TO  AVOID  DIVISION  BY  ZERO" 

L2=L 1-1 


IF  L 2 <  =  L 2 m ax  THEN  518 

PRINT  " L2  REDUCED  TO  L2max,  TO  AVOID  WRAP-AROUND  REGION 
L  2  =  L  2  m  ax 

IF  N2>L2+L2  THEN  558 

PRINT  " N 2  INCREASED,  SO  AS  TO  BE  GREATER  THAN  L2*2n 

N 2  m  = I  NT (LOG ( L2  +  L2+ 1 ) /LOG  < 2 ) >  + 1 

N2=2AN2rn 

PRINT  ,!L2  = "  ;  L2 ;  "  N2  =";N2 

N  2  rn  =  N  2  - 1 

RED  I  M  X  <  0 :  N2rn  )  ,  Y  <  0  :  N2rn  > 
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580 

A  =  P I  /  L  2 

590 

X ( 0 ) =Xa<  0  > 

60O 

FOR  N s  =  1  TO  L2 

610 

B=Ns*A 

620 

W=fl0+Rl*COS<B>+A2*COS< 

630 

X < Ns > =  X  <  N2-Ns > =  Xa< Ns ) * 

648 

NEXT  Ns 

650. 

FOR  N s  =  L 2  + 1  TO  N2m-L2 

660 

X  <  Ns ) =0 . 

670 

NEXT  Ns 

6S0 

MAT  Y  =  < 0 . > 

690 

CALL  Fftl4z<N2,X<*>,Y< 

700 

MAT  X  =  X  /  (  F  L  T 00*L1*N1* 

710 

Bi g=MAX<X<*> ) 

720 

Db=LGT<Big>*10. 

730 

PRINT  “Big  =  " ; Bi q;  “ 

740 

PLOTTER  IS  "GRAPHICS" 

750 

GRAPHICS  ON 

760 

WINDOW  0. , N2/2, -70. , 0. 

770 

GRID  N2/3, 10. 

780 

FOR  N s = 0  TO  N2/2 

790 

W  =  M  A  X  <  X  < N  s ) ,  1 . E-20  > 

800 

PLOT  Ns, L G T <  W  )  *  1 0  • - D b 

810 

NEXT  Ns 

828 

PENUP 

330 

GOTO  430 

340 

END 

350 

! 

360 

SUB  Fftl4z< DOUBLE  N , R E 

8  7  0 

!  This  SUB  is  listed  ab o 

1670 

SUBEND 

1680 

! 

1690 

SUB  Data< DOUBLE  Ks,Ll, 

1700 

DOUBLE  J 1 , J 2 , J 

1710 

W 1 =2 . *P I  * . 25 

1720 

W  2  =  2 . *PI*. 265 

1730 

T  1  =  .  3  1 

1740 

T2=. 77 

1750 

D  b  2 = -  3  0  . 

1760 

A2= 1 0 . A<Db2/20. > 

1770 

J2=Ks*Ll-l 

1730 

J1=J2-L1+1 

1790 

FOR  J  =  J 1  TO  J2 

1300 

P1-W1*J+T1 

1310 

P2=W2*J+T2 

1820 

X< J-Jl  )=C0S<P1 >+A2*C0S 

1830 

!  RND-.5  has  power  1/12 

1340 

NEXT  J 

1350 

SUBEND 

LfiG  WEIGHTING 


INTO  FREQUENCY  DOMAIN 


N2> 


dB  max  =  ";Db 


N  =  2-^1  N T E G E R <  = 2  M  4  =  1  6 3 8 4  ;  0 


REAL 


;  /  *  )  ) i 


TONE 

TONE 

TONE 

TONE 


FREQUENCY 
FREQUENCY 
PHASE  1 
PHASE  2 


RELATIVE  TONE  STRENGTH 


■<RND-.5>*.316 
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